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Bilayer graphene has been predicted to give unprecedented tunability of the electron-electron 
interaction with the help of external parameters, allowing one to stabilize different fractional quan¬ 
tum Hall states. Recent experimental works make theoretical analysis of such systems extremely 
relevant. In this paper we describe a methodology for investigating the possibility of realizing spe¬ 
cific fractional quantum Hall states in bilayer graphene taking into account polarization effects and 
virtual interband transitions. We apply this methodology to explore the possibility of realizing the 
Moore-Read Pfaffian state in bilayer graphene. 
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Introduction 


I. NUMERICAL DIAGONALIZATION 
APPROACH TO THE QUANTUM HALL 
EFFECT IN NON-RELATIVISTIC SYSTEMS 
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Since the discovery of graphene [T] it has been well- 
understood that this material opens new horizons on the 
investigation of quantum Hall effects (QHE). The reasons 
for this are a much better confinement of the electron 
gas in a 2D plane, longer mean free paths, and larger 
cyclotron gaps than can be achieved in GaAs quantum 
wells. The integer and fractional QHE have been ob¬ 
served in systems with suspended graphene and graphene 
on different substrates HU. The sequence of conduct¬ 
ing plateaux in graphene was found to differ from that in 
GaAs, which is due to a peculiar structure of the single 
particle spectrum. 

Bilayer graphene (BLG) has attracted attention of the 
fractional quantum Hall effect (EQHE) community for 
some time. The main reason for that is that this mate¬ 
rial allows for unprecedented tunability of the parameters 
important for EQHE by means of external electric and 
magnetic fields mm- It turns out that the effects of Lan¬ 
dau level mixing, vacuum polarization etc. are extremely 
important in this system [7], making it complicated to 
analyze theoretically. However, recent experimental ad¬ 
vances in observation of EQHE in BLG [SUTO] call for 
theoretical analysis. This paper expounds in detail a mi¬ 
croscopic study of the possibility of realizing the so-called 
Moore-Read Pfaffian EQHE state m in BLG, which was 
briefly reported in [ 7 ]. 

We perform the analysis of Landau level mixing, vac¬ 
uum polarization and some other effects in order to es¬ 
tablish the conditions under which theoretical predictions 
are reliable. Under appropriate conditions, we treat level 
mixing pertubatively in the spirit of nans] and inves¬ 
tigate the possibility of realizing the Moore-Read state 
using exact numerical diagonalization. The methodology 
presented here can be readily applied to study the pos¬ 
sibility of realizing any other EQHE state in BLG given 
the state’s trial wave function. 


In this section we recall how Landau levels (LLs) 
emerge in a two-dimensional non-relativistic system in 
a magnetic field. We also introduce some notation that 
will be used in the following sections. 


I.l. Problem of a single non-relativistic electron in 
a magnetic field 


The single electron Hamiltonian in the uniform mag¬ 
netic field, perpendicular to the plane of the system, is 

TT^ 

^l-part = ( 1 ) 

where tt = tt^) (since the system is two-dimensional), 
TTi = Pi + eAijc^ Pi = —ihdi^ e is the elementary charge, 
A = [B X r]/2 is the vector potential of the uniform 
magnetic field B = —Bbz, Sz is the ^-component of the 
electron spin, rUe is the free electron mass, m* is the 
effective mass of an electron, g is the Lande p-factor, and 
Pb = e/i/(2mec) is the Bohr magneton. 

We introduce the magnetic length /, the cyclotron fre¬ 
quency cjc, and the complex coordinate w in the plane: 


7 / 


eB 


X + iy _ X — iy 
w = —^—. 


m*c ’ I l 

We also introduce operators d, d'*', 6, 

a = V2(a+|), at = V2(-5+i), 
6 = V2(a+i), 6t = V2(-a+|). 


( 2 ) 

( 3 ) 

( 4 ) 


where d and B denote d/dw and d/dw respectively. All 
commutation relations between these four operators are 
trivial except for the following: 


The first section reminds the reader of the method- 
oiogy for studying EQHE in conventionai (”non- 
reiativistic”) systems based on exact numericai diagonai- 
ization of the system Hamiltonian in the single Landau 
level approximation. 

In the second section we describe a formalism which 
allows the use of the same methodology in bilayer 
graphene. 

In the third section we discuss the effects of inter¬ 
action of different Landau levels and incorporation of 
them into the single-Landau-level-based methodology in¬ 
troduced earlier. 

Einally, in the fourth section we apply our methodol¬ 
ogy that takes into account the deviations from the single 
Landau level approximation to study the possibility of 
realizing the Moore-Read m state in bilayer graphene. 


[a, 



( 5 ) 


We can then rewrite the Hamiltonian in the form 




1—part 


= huOr 


a^a ■ 


- ] - gpBBSz^ 


(6) 


The operators d, are similar to the ladder operators 
in the problem of the harmonic oscillator, thus the sys¬ 
tem’s spectrum consists of Landau levels with energies 
En = hLJc{n + 1/2) — gpBBSz^ n G Sz = ±1/2. Op¬ 
erators 6, w commute with the Hamiltonian, thus they 
transform one state into another state within the same 
Landau level. 

Let us consider the operator of the 2 :-projection of the 
orbital angular momentum: 

L = Lz/h = [r X p]^ jh = zd — zB = — d)d. (7) 
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One can see that 


L, a 

= 

L,a1' 

= 

L,b 

= -b, 

L,b^ 


-^5 —part 


= 0 . 


( 9 ) 


Thus the eigenstates of the Hamiltonian 0 can be 
labeled by three quantum numbers: the Landau level 
number n G the 2 :-projection of the angular momen¬ 
tum m G (Z+ — n) and the electron spin projection 
= ± 1/2 


\n,m,Sz) = tpnmiw) 'S’lSz). (10) 

The orbital wave function ^pnm can be expressed as 

1 


'lpnm{w) = 


ipooiw) = 


■\/n!(n + m)\ 


( 11 ) 


1 _ 

e 4 ^ a'tpoo = Hoo = 0. (12) 




Note that the wave functions (11) are polynomials of 
complex coordinates multiplied by the exponential 

exp (—|rep/4) which is the same for all of the states. In 
particular, 


'0n=O,m ('^) 


1 



4 


(13) 


So, the system’s energy levels are Landau levels with 
energies — ^c(^±l/2)—with eigenstates 

in a LL labeled by the angular momentum projection 
m > —n. 


In a finite sample there is only a finite number of states 
available to an electron in a Landau level. One can esti¬ 
mate their number using the fact that the states 
are spatially localized: the number of states in a Landau 
level of a finite round sample is approximately equal to 
the number of states (11)-([T^ that are localized mainly 


in the area of the sample. Thus, one can introduce the 
filling factor u: 


V = 2TTfne = Ne/Norh, (14) 

where Ue is the density of electrons, is the total num¬ 
ber of electrons in the system and A^orb is the number of 
orbitals in a LL in the sample. 

Typically in GaAs heterostructures m* ~ O.OTme, 
g ^ —0.4, thus LLs with the same number n but different 
spin projections form closely spaced doublets. For a typ¬ 
ical fractional quantum Hall (FQH) experiment in such 
systems the characteristic Coulomb energy scale /{nl) 
{n is the dielectric constant, in GaAs n 13) is on the 
same order as the cyclotron frequency hwc. For H > 5T 
the interaction energy scale is less than the inter-doublet 
spacing. Therefore, in these conditions the Coulomb in¬ 
teraction cannot throw the electrons from a Landau level 
to other Landau levels efficiently. 

In this case for any filling factor one can expect only 
one doublet to be partially filled with others either fully 


filled or completely empty. Thence, it is not too bad an 
approximation to restrict consideration to the electrons 
in the partially filled doublet. Corrections to this pic¬ 
ture can be taken into account by means of perturbation 
theory [MIT]. However, in this section we neglect them. 

In a partially filled doublet, the lowest energy is usu¬ 
ally achieved when the electrons form a spin-polarized 
state. In such a state only one of the two Landau lev¬ 
els in a doublet is partially occupied, and the other is 
either empty or fully occupied. There are two reasons 
for that. One is that such states minimize the Coulomb 
interaction exchange energy (if the interaction potential 
decreases monotonically with distance, which is typically 
the case). Another reason for the electrons to form spin- 
polarized states is the Zeeman splitting (even though it 
is small). 

In the remaining part of this section we only consider 
one Landau level, neglecting the influence of other Lan¬ 
dau levels. This approximation is called the single Lan¬ 
dau level approximation (SLLA). We also assume that 
the state is spin-polarized, therefore we suppress the spin 
variables. 


1.2. Interaction of two electrons in a 
non-relativistic Landau level 

We begin the discussion of the many-body problem 
with the two-particle case. For interaction potentials 
which depend only on the distance between the electrons 
this problem can be solved exactly. This solution gives 
an opportunity to introduce some important notions. 

The two-electron Hamiltonian can be written as fol¬ 
lows: 

H2 —part Hhee + y{r), (15) 

l^free “ — part,l — part, 2 5 (16) 

where r = |ri — r 2 | = l\wi — 1 ^ 2 !, V{r) is the interaction 
potential, e.g., the Coulomb potential. 

Since we are working in the SLLA approximation, the 
single-particle part of the Hamiltonian is proportional 
to the identity operator and can be excluded from the 
consideration. Thus to diagonalize the Hamiltonian we 
only need to diagonalize the interaction potential opera¬ 
tor V (r) in the Hilbert space spanned by vectors 

|mi,m2) = ® 1 ^ 2 ) - |m2) 0 |mi)), (17) 

with the angular momenta of the two electrons, mi and 
m 2 , taking all the possible values in the LL considered. 

We introduce ^-projections of the relative angular mo¬ 
mentum and the angular momentum of the center of 
mass: 


^rei = ( ^[(ri -r2) X (pi 


■P2) 


— ^±1 -j- L 2 — b\b2 — ^ 2^1 ± ci\ci2 ± , (18) 
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Lem = ( ^[(ri + ^2) X (pi + P2) 

— + L2 + b\b2 + ^ 2^1 ~ CL1CL2 ~ cl 2 ^i^ • ( 19 ) 

These operators projected onto a single LL have the 
following form: 


iu = l( 


— 7: {Li ^ L2 — b 


I &2 - blh) 


Llm 


- ^Li + L 2 + b\b2 + &2^i) • 


( 20 ) 

( 21 ) 


The raising and lowering operators for this ’’single-level 
angular momenta” are b\ ^ b^ and hi ^ 62 respectively. 
With the help of these operators we can represent the 
eigenstates of the ’’single-level angular momenta” in the 
(n^Sz) Landau level as follows: 


|m, M) = 


1 


press Vm 

\/2'"+^777!M! 


form [B?]: 

(6l - blribl + bl)^ii^n,-n)l{i’n,-n)2, 

(22) 



V{q) = 

1 

II 

(23) 


LUm,M) = iM-n)\m,M). 

(24) 



Here M, m > 0. We will say that |m, M) is a state with 
the relative angular momentum m and the center of mass 
angular momentum M. Since every state has to be anti¬ 
symmetric under the permutation of the electrons, only 
states with odd m are present in our Hilbert space. The 
states |m,M) for m G 2Z+ + 1 and M G form a 
complete orthonormal basis. 

Commutation relations of the ” angular momenta” with 
the operator V = V{r) are 




rel/cm 


= 0 , 


V,b\ + bl 

= 0 , 

V,bi+b2 

'V,bl-bi 

^ 0 , 

'V,h-b2 


0 , 


(25) 

(26) 
(27) 


In the following section we shall also need a more gen¬ 
eral matrix element 

^ - x 

' V 2 -+^m!M! 

{b[ - bl)^{b\ + bl)^{ipn,-m)ii'ipn2- 1 x 2 ) 2 , (29) 

= ((m, 712, m, My Im, 712, m,M)) = 

((771,772, 777, 0y|77i, 772, 777,0)). (30) 

It is easy to check that for the po- 

tentials depending on the distance between the electrons 
only. Note that for technical reasons we do not impose 
antisymmetry on the wave functions ( [ 2 ^ , to emphasize 
this we mark these states with a double bracket )). 

For computations, it is often more convenient to ex- 




27r 

7^ 


Yirii,n2) 




d^rV{r)e-^^^l^ = 

POO 

/ V{r)Jo{qr/l)rdr = 

Jo 

pOO 

27r / V{lx)Jo{qx)xdx, (31) 
Jo 


Viq)Lmiq^)L„, {q^/2)L^, (gV2)e-^' 

(32) 

where Jq is the zeroth Bessel function of the first kind, 
are the Laguerre polynomials, I is the magnetic length. 
Derivation of this formula is presented in Appendix 
Thus, the electrostatic interaction between the elec¬ 
trons located in one Landau level can be expressed 
through a countable set of pseudopotentials , where 

n is the LL number, and m G 2Z+ + 1 is the ’’relative 
angular momentum” of the two interacting electrons. 


1.3. Many-particle problem 


Thus the interaction potential operator can be repre¬ 
sented in the LL as 

V{r) = Y, \m,M)V^’^Hm,Ml (28) 

m,M 

which solves the two-body problem in the LL. 

Matrix elements which parametrize the opera¬ 

tor are called pseudopotential coefficients (or just pseu¬ 
dopotentials), they were first introduced in [18]. The 
connection of the pseudopotentials with the interaction 
potential’s matrix elements is obvious since the states 
|m, M) are orthonormal. 


Here we discuss the problem of many electrons in 
a Landau level and the numerical diagonalization ap¬ 
proach. 

Since we know how to express the Hamiltonian of the 
electron-electron interaction in the Hilbert space of two 
electrons in a Landau level, we can, in principle, ex¬ 
press the many-particle system’s Hamiltonian through 
the pseudopotentials. Typically, such a Hamiltonian can¬ 
not be solved analytically, therefore numerical diagonal¬ 
ization is used. A standard complication is that the 
Hamiltonian is an infinite matrix (since there are an infi¬ 
nite number of orbitals in a Landau level), while numer¬ 
ical diagonalization can only be used for finite matrices. 
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Therefore, several strategies are used to restrict the 
system size. One is to restrict the maximum orbital num¬ 
ber available to the electrons (i.e. consider only those 
states in the LL where electrons occupy orbitals with an¬ 
gular momentum quantum number m < rUmax)- Such a 
model is called the ’’hard cutoff’ model for a system on 
disk. It is also possible to restrict the total angular mo¬ 
mentum of the electrons = Mmax), which gives 

the ’’soft cutoff’ model for a system on disk [2QH^ . The 
third widely used method is to consider a ’’system on 
sphere” [18] (two-dimensional finite sphere with the uni¬ 
form magnetic field transverse to the sphere is considered 
instead of plane). A Landau level in a system on sphere 
is finite (has a finite number of orbitals) from the very 
beginning, so one does not need to introduce an artifi¬ 
cial boundary. The wave functions of the single particle 
states and pseudopotentials are expressed in a somewhat 
different way (so the matrix of the Hamiltonian is ex¬ 
pressed somewhat differently via spherical pseudopoten¬ 
tials). However, for large enough systems the results on 
sphere should coincide with the results on disk (since the 
curvature of the sphere plays little role then). That’s 
why planar pseudopotentials are often used for diagonal- 
ization on sphere (see e.g. [23]). In this work we do 
a similar thing: we use diagonalization on sphere with 
planar pseudopotentials. 

The Hamiltonian of a system on sphere/disk is a finite 
matrix that can be expressed in terms of pseudopotentials 
introduced in the previous subsection and diagonalized 
numerically. This enables us to find the spectrum and 
the eigenstates of the system. A typical thing to do then 
is to compare the numerically found ground state (and, 
possibly, the excited states) with some trial wave function 
to check whether the real state is close to the proposed 
trial state.^ 

There is some peculiarity in choosing the number of 
electrons and orbitals in the system. If one studies the 
filling factor z/, then by definition in the thermodynamic 
limit number of electrons in the LL considered and the 
number of orbitals available to them A'orb are related by 
A^e/Aorb ~ where {u} denotes the fractional part of 
the filling factor. On the contrary, trial wave functions 
(as can be seen from examples below) fix the relation 
between the two numbers not approximately but exactly: 


properties of the system in the thermodynamic limit to 
be independent of the precise ratio between the num¬ 
ber of electrons and the number of orbitals; but in order 
to compare an exact state with a trial wave function, the 
number of orbitals and electrons in each of the two states 
should be related as described by Eq. (33). 

So, the procedure of numerical finding the system’s 
ground state and its comparison with trial state is as 
follows: choose the trial state with which to compare; 
choose the number of electrons and orbitals in such a way 
that it corresponds to the trial state; find pseudopoten¬ 
tials; calculate the system’s Hamiltonian and diagonalize 
it; calculate the scalar product of the numerically found 
ground state with the trial state (the closer it is to 1 the 
more similar the states are). 

Usually the numerical diagonalization can be per¬ 
formed only for relatively small numbers of electrons 
(around 10 to 20) in most cases. This is far from the ther¬ 
modynamic limit. However, it has historically tended to 
be the case that even systems this small can be fairly rep¬ 
resentative of the thermodynamic limit [TO] [24] . Strictly 
speaking, one should attempt extrapolation to infinite 
size. However, this is beyond the scope of our project 
and we do not believe that our main result would be 
qualitatively changed. 

Before proceeding to application of this method to bi¬ 
layer graphene, we show several examples of trial wave 
functions in the next sub-subsection. 


L3a. Examples of trial wave funetions 

Here we consider several examples of trial wave func¬ 
tions in order to understand how they look and how to 
interpret them (for the purposes of numerical diagonal¬ 
ization). For simplicity, we present the trial wave func¬ 
tions for a system on disk. 

The simplest example is a trial wave function for the 
fully occupied zeroth LL. Let N be the number of elec¬ 
trons which occupy the first N orbitals of the n = 0 level. 
Due to the Pauli principle the only possible state is the 
Slater determinant of all the occupied single-particle or¬ 
bitals: 


A^orb = A^e/W-5 + l. (33) 

Number S is called ’’shift” and can be different for dif¬ 
ferent trial wave functions.^. Of course, we expect the 


-IpilUi, ...,Wn) 


det i’n=0,m{u>i) OC 

l<i<Ar,0<m<JV-l 


(X 


1 Wi 
1 W2 


N-1 

1 

N-1 


(34) 


^ There is a correspondence between trial states which are pro¬ 
posed for the sphere and for the plane, so a result of the diago- 
nalization on a sphere can be compared with a trial state just in 
the same way as a result of the diagonlization on a disk. 

^ The summand +1 is for the number of the last available orbital 
in the zeroth LL mmax in a system on disk to be expressed as 
nimax = Ne/ {v'} — S. This is a commonly used definition of the 
shift. 


1 Wn • • • w 


N-1 

N 


where we used the explicit form of the single-particle 
wave functions in the zeroth LL (13). The determinant 


on the r.h.s. is the well known Vandermonde determi¬ 
nant. Therefore, we can write down the answer for the 
wave function, which, up to the normalization constant 







6 


A/", looks as follows: 

'iP(wi,...,wn) JJ {Wi-Wj) = 

l<i<j<N 

(35) 

This example illustrates the fact that any wave func¬ 
tion of electrons in the zeroth LL can be expressed as a 
polynomial of coordinates Wi — no Wi — times the Gaus¬ 
sian weight. Below in this sub-subsection instead of the 
wave function ...^wn) we will write out the polyno¬ 

mial part P{wi, ...,wn)- For example, for the fully filled 
zeroth LL 

P{wi,...,wn) = (Wi-Wj). (36) 

If electrons don’t fill the whole LL they will try to keep 
as big a distance from each other as possible (because of 
the Coulomb repulsion). Starting from this argument, 
R. Laughlin proposed his famous trial wave function for 
the filling factor z/ = 1 /3 [25] : 

P{wi,...,wn) = (Wi-Wj)^. (37) 

l<i<j<N 

It is easy to convince oneself that it indeed corresponds 
to z/ = 1/3 by counting the number of orbitals used by 
the electrons in this wave function. It has the shift S = 3 
(in contrast to the full filling, where S = 1). The key idea 
is that the power 3 significantly reduces the probability 
of finding two electrons close to each other. 

This wave function has been generalized for the fillings 
u = 1/m: 

P{wi,...,wn) = tJ {Wi-Wj)'^. (38) 

However, since the wave function of the electrons should 
be antisymmetric, m has to be odd. So, this wave func¬ 
tion can be used only for fillings with odd denominators. 

There had not been any need in description of even de¬ 
nominators until the z/ = 5/2 FQHE was observed [26] . 
Moore and Read in 1991 proposed their trial wave func¬ 
tion for a half-filled LL^ m- The wave function, if writ¬ 
ten for the zeroth LL, looks like 

P(n;i, ...,n;Ar) = Pfaff {wi-Wjf = 

\W^ Wj J ;L<i<j<Ar 

. .n ^ 1 1 1 A 

AntiSymm-...- x 

\Wi -W2Ws- W4 Wn-1 - wn J 


^ As it has been mentioned already, it is assumed that only the 
partially filled level is important. 


AntiSymm (...) denotes the antisymmetrization opera¬ 
tor. The antisymmetrized combination which is present 
here is the Pfafhan of the matrix Mij {Mu = 0, Mij = 
1 /{wi — Wj)). After this expression the Moore-Read wave 
function is also often called the Pfafhan state. The hll- 
ing factor associated with this wave function is u = 1/2, 
and the shift is S' = 3. The wave function is evidently 
antisymmetric. 

One can write wave functions for higher Landau levels 
in a similar explicit fashion, but, in fact, for numerical 
comparison one only needs the coefficients of the state 
vector expanded in the basis of orbital occupation num¬ 
bers. For a zeroth LL wave function those coefficients 
can be obtained from the polynomial representing it by 
expanding the polynomial into a linear combination of 
monomials — each wf up to the normalization factor 
corresponds to an electron occupying the state 7 /^= 0 ,m=/c- 
One can also interpret the zeroth LL trial wave function 
as a wave function for a higher LL. For that one should 
replace 7 / 0 ,/c ^ '^n,k-n f he very end of the procedure of 
getting the coefficients. Therefore, polynomials of elec¬ 
trons’ coordinates Wi are used for representing trial wave 
functions for both zeroth LL and the higher LLs. 

Similarly, one can map angular momentum orbitals 
of the zeroth LL of a planar non-relativistic system to 
the corresponding angular momentum orbitals of a LL 
of BLG (or of the system on sphere), so that any anti¬ 
symmetric polynomial (zeroth LL wave function) can be 
used as a trial wave function of any LL in BLG or the 
system on sphere. 

Thus, all the trial states, including the Moore-Read 
Pfaffian, can be used for higher LLs. In the n = 1 LL the 
Moore-Read state’s overlap with the numerically found 
ground state for 12 electrons is close to 0.7^. This is not 
as impressive as Laughlin’s 98 — 99%, but still remarkably 
good for the Hilbert space of dimension over 16 thousand 
(two random vectors would have an overlap near 1/16000 
in such space). The FQHE with u = b/2 (which corre¬ 
sponds to a half-filled n = 1 LL) is observed in GaAs 
heterostructures. Numerical studies strongly sug¬ 

gest that the state is either the Moore-Read Pfaffian or 
its particle-hole conjugate called anti-Pfaffian [32l l33] . 
Recent experimental studies seem to rule out the Pfaf¬ 
fian |34l|35], however, whether the state is anti-Pfaffian 
remains to be confirmed. One difficulty in dealing with 
the z/ = 5/2 fraction in GaAs is the extreme fragility of 
the corresponding state. 


1.4. Summary of the section 

In this section we review the basis of numerical diag- 
onalization methodology for non-relativistic EQHE sys¬ 
tems: introduce Landau levels, briefly discuss the ap- 


^ By overlap we mean scalar product’s absolute value squared. 
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plicability of the SLLA to the GaAs heterostructures, 
introduce pseudopotentials, and discuss peculiarities of 
numerical diagonalization in the non-relativistic systems. 
We also discuss several examples of trial wave functions, 
including the Moore-Read Pfaffian, and their representa¬ 
tion in the form of holomorphic polynomials of complex 
coordinates. 


II. NUMERICAL DIAGONALIZATION 
APPROACH TO THE QUANTUM HALL 
EFFECT IN BILAYER GRAPHENE (SINGLE 
LANDAU LEVEL APPROXIMATION) 


In this section we discuss peculiarities of the numer¬ 
ical diagonalization method in bilayer graphene in the 
SLLA. Landau levels in bilayer graphene are introduced, 
expressions for the pseudopotentials are derived. 


H.l. Bilayer graphene. Hamiltonian of a free 
electron in bilayer graphene 


Graphene is a one-atom thick layer of graphite, or in 
other words — a two dimensional honeycomb lattice of 
carbon atoms. Bilayer graphene (BLG) is formed of 
two layers of graphene (two graphene sheets) with cer¬ 
tain matching of lattice points. For a detailed review 
on graphene and bilayer graphene see Ref. [36]. We are 
going to recall only the facts necessary for the following 
consideration. 

The Fermi surface of undoped graphene consists of two 
points (valleys, usually denoted as K and K') in the first 
Brillouin zone. One can usually neglect jumping of elec¬ 
trons between the valleys.^ Therefore, the valley index 
of an electron is a good quantum number. 

In BLG the low-energy excitations are also located 
near the same Fermi points in the momentum space. The 
low-energy BLG Hamiltonian (without external magnetic 


field) can then be written as [36ll37] ^ 
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(40) 


where ^ = ±1 is the valley index such that ^ = +1 corre¬ 
sponds to K and ^ = — 1 corresponds to A', tt = Px+iPy 
is the complex momentum. The spectrum has a mini-gap 
21 /, which can be tuned by the external electric field per¬ 
pendicular to the bilayer graphene sheet.^ We will call U 
the ’’mini-gap parameter”. The Fermi velocity is taken 
to be u ^ 10 ^ m/s, and the interlayer hopping constant 
is taken to be 71 « 0.35 eV [36] . 

In the absence of the external electric field (when 
U = 0) the low-energy spectrum has quadratic form E = 
±|7rp/(2m*), where the effective mass m* = 71/(2u^) 

0.03me [371 . 


H.2. Problem of a single BLG electron in a 
magnetic field 

The Hamiltonian of an electron in bilayer graphene in 
the perpendicular magnetic field is obtained by making 
the derivatives covariant and taking the spin energy into 
account: 


^BLG _ 
1—part 
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utt’*' \ 
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U 

UTT 
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utt’*' 
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y UTT 
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Gi 
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(41) 


where ^ = ±1 is for the two valleys, tt = tTx + i'Ky (see 
definition of TTj after formula Q). In BLG the Lande 
factor g ^ 2. Without loss of generality we will consider 
only the case B^U > 0 . 

It is easy to express the complex momenta through the 
operators p||4| ): 


TT = —iV2hl ^d, TT^ = iV2hl . (42) 


H 


Thus the Hamiltonian can be expressed as 
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i'ya) \ 
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u 

—i^a 

u 
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\- 

-iyd 

0 


-u ) 


- qPbBSz, 


(43) 


^ This is due to the fact that jumping needs transfer of a quite big 
momentum (of the order of h/a, where a is the lattice constant 
and has value around 0.25 nm). For example, matrix element 
of the Coulomb potential decreases like 1/q a,s the transferred 
momentum q grows. Therefore, jumping between the valleys is 
suppressed, with controlling parameter being the ratio of the 
lattice constant a to the typical spatial scale one is interested in. 
(In our case this is the magnetic length 1; typical values of the 
magnetic length are around I = 10 nm). 


® The Hamiltonian is written in the basis corresponding to the 
atomic sites A, B, A, B in the K valley and JB, A, B, A in the 
K' valley. The sites A and B are situated in the bottom graphene 
layer, while A and B are in the top layer. Our convention is the 
same as the one used in Refs. [selllT] except for a redefinition of 
U. 

^ One can think that the electrostatic potential of one layer is B, 
while the other layer’s potential is —U. 










where we introduced ujc = eBj^rn^c) = 2v^eBj^\c (after 
definition of ref. ISZ1)> 7^ = 7i/(^c), W = Ul{hu}c)- 
This Hamiltonian does not commute with the 2 ;- 
projection of the orbital angular momentum L defined 
in Eq. 0 . We introduce the ^-projection of the ’’pseu¬ 
dospin angular momentum” E: 


/I 

0 
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0 

-1 
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0 
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Vo 
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0 

0/ 


(44) 


Then the 2 )-projection of the full angular momentum J = 
L + E does commute with the Hamiltonian. 


Now it is easy to express the general form of the spatial 
part of the Hamiltonian’s (43) eigenstates through the 
non-relativistic wave functions ( TT][T2 ): 


^ nm 


( \ 

2,m+2 1 

l,m+l I 

\ — J 


(45) 


the sense of the number n is similar to the Landau level 
number, while m corresponds to the projection of the 
full angular momentum 2z — rn -\- 1. The amplitudes 
Bn^ Cn, do uot depend on m. 

Acting on this wave function with the Hamiltonian and 
demanding it to be an eigenfunction we find the equation 
for the eigenvalues: 

{{u-Cef—y‘^n){{u+^ef—j‘^(n-l)) (46) 


where e = {E ^ 2fiBBSz)/{huJc)^ and E is the energy. 

Finding the single particle spectrum for the realistic 
values of the parameters, we see that the levels split into 
two groups: the one with \E\ < 71 and the one with 
\E\ > 7 i. The Zeeman splitting is negligibly small, just 
like the non-relativistic case. The levels with \E\ < 71 , 
which we are interested in can be characterized by five 
quantum numbers: the valley index the LL number 
n G the full angular momentum projection = m+1 
with m G (Z+ — n), s = ±1 (which shows whether the 
energy is positive or negative) and Sz = ±1/2. Thus the 
wave functions (their spatial components) in the n-th LL 
look like 




/ \ 

V^n—l,m+l 

\ l,m+l / 


(47) 


The amplitudes which are present in this formula can be 


expressed as follows 


^ _ Vn-l {u - ^e «^)2 - 


u + ^e' 




= - 


v 2^ 


V, 


^73 


V, 


Dis ^ 

70^ 


(48) 

(49) 

(50) 

(51) 


where V is a normalization constant. Obviously, they de¬ 
pend on the magnetic field B and the mini-gap parameter 

U. 

Before considering the two particle problem in bilayer 
graphene, we have a look at the single-particle spectrum. 
Fig. shows the dependence of the several lowest Lan¬ 
dau levels on the magnetic field for U = bO meV. Only 
the positive part of the spectrum is shown, the negative 
part can be obtained with the help of electron-hole con¬ 
jugation = — 5 “^’^). Each positive LL is labeled 

by a pair of quantum numbers (n,^). When we need to 
work with negative LLs, we label them with j = (n, s). 

One can see that for large values of the magnetic field 
the levels form quasidegenerate doublets which are sep¬ 
arated by energies of the order of huOc- Fig. shows 
the dependence of the same LLs’ energies on the mini¬ 
gap parameter U for the magnetic field H = 10 T. Note 
that for large enough values of U (or for small enough 
values of B) multiple crossings of Landau levels occur. 
It is easy to understand that when several LLs are close 
to each other (for small magnetic fields/large mini-gaps 
when the LLs cross, or for small values of the mini-gap 
when the levels in a doublet are almost degenerate) sig¬ 
nificant deviation from the SLLA can occur. Thus, the 
applicability of the SLLA puts constraints onto the ex¬ 
ternal parameters. This point is discussed in details in 
section Ell 


II.3. Interaction of two electrons in a Landau level 
of bilayer graphene 

The two-particle problem within the SLLA in BLG can 
be solved analogously to the non-relativistic case. 

Define the projections of the relative and the center of 
mass full angular momenta to a single Landau level: 

±i/cm = ±i/cm + 2(^1 + ^2) = 

— (Ji -|- J 2 T (^1^2 + b^bi)). (52) 

The commutation relations of the projected angular 
momenta and their raising and lowering operators with 
the parts of the two-particle Hamiltonian are similar to 
the non-relativistic case: 

= H^h^ + V{r), (53) 

^BLG ^BLG I ^BLG fKA\ 

-^free -^1 —part,l ' -^1 —part,2? 
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Figure 1: Dependence of the lowest LLs’ energies on 
(a) the magnetic field B for = 50 meV, and (b) on 
the mini-gap parameter U for B = 10 T, Each level is 
labeled by a pair of quantum numbers (n, ^). Only positive- 
energy part of the spectrum is shown. (Color online). 


^BLG jP 
^free 5 '^rel/cm 

= 0 , 


(55) 


= 0 , 

V 

’ ^rel/cm 

= 0, 

(56) 


5> 

+ 

to—H 

= 0 , 

vM + k 

= 0, 

(57) 


'v,b\-bi 

7^ 0 , 

'vM-k 

^0. 

(58) 


The simultaneous eigenstates of the two ” angular mo¬ 
menta” (52) in a Landau level have the form 


|m,M) 


1 

“ V 2 ™+^m!M! 

{b\ - biribi+bi)^{^i:-nh{^i:-n)2, 


(59) 


M) = {m - n + l)\m, M), (60) 

JPJm,M) = {M -n+l)\m,M). (61) 

Thus, just like in a non-relativistic system, the two- 
particle interaction potential in a Landau level can be 
expressed through pseudopotentials: 

V{r) = \m,M)VZ’^’^m,M\. (62) 

m,M 


The expression for the pseudopotentials in terms of ma¬ 
trix elements of the potential is straightforward as the 
states |m, M) are orthonormal. These pseudopotentials 
can be expressed via the non-relativistic pseudopoten¬ 
tials: 


Yn,^,s 
^ m 




(63) 


This expression makes obvious the possibility of tuning 
of the pseudopotentials by changing the values of the am¬ 
plitudes Aff ^ ^ . Since the amplitudes depend 

on the external parameters U and 5, the pseudopoten¬ 
tials can be tuned with the help of external perpendicular 
electric and magnetic fields. 

For practical purposes it is useful to incorporate for¬ 
mula (32) into Eq. (1^ which leads to 


Yn,^,s _ 


V{q)Lra{q^) (Fy(gV2))'e-^'^, (64) 


Fy(gV2) = |4l'in(9V2) + l4l'in-2(gV2) + 

(|4^|^ + |4*|^)i:„_i(gV2). (65) 

After pseudopotentials are calculated the numeri¬ 
cal diagonalization procedure is identical to the non- 
relativistic case^. 

We emphasize that the analysis in this section is done 
strictly within the SLLA. As is discussed in the next sec¬ 
tion, in BLG the SLLA is less justified than in GaAs 
systems. There are, however, conditions under which the 
SLLA is a good approximation. In the latter case the 
effects of other LLs can be incorporated into the SLLA 
by means of perturbation theory. 


II.4. Summary of the section 

In this section the explicit formulae for the SLLA in 
BLG are provided (wave functions in a Landau level, ex¬ 
pression for the pseudopotentials). It is shown that appli¬ 
cation of the numerical diagonalization methodology to 
BLG system within the SLLA does not differ too much 
from the application to a non-relativistic system. 


Recall that though the trial wave functions are written in the 
form of complex polynomials, they, in fact, give decomposition 
of the wave function into Slater determinants of the single parti¬ 
cle states. Thus they are applicable to any system with Landau 
levels having structure similar to the non-relativistic case, so they 
are applicable to the BLG. If a trial state is written in the basis 
of occupation numbers the only difference to the diagonaliza¬ 
tion and comparison procedure is that one has to use the BLG 
pseudopotentials instead of the non-relativistic ones. 
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III. DEVIATIONS FROM THE SINGLE 
LANDAU LEVEL APPROXIMATION IN 
BILAYER GRAPHENE 


As discussed in the previous section, studying the 
FQHE in BLG within the SLLA does not differ too much 
from the non-relativistic case. However, as it has already 
been mentioned the single-particle spectrum of BLG re¬ 
stricts the applicability of the SLLA, and under the usual 
experimental conditions the restriction is significant. De¬ 
pending on several factors, there are three regimes: 

1. the SLLA is fully reliable; 

2 . the SLLA is completely inapplicable because one 
has to consider several Landau levels together since 
the electrons partially occupy each of those; 

3. intermediate regime, the effects of the presence of 
other Landau levels can be incorporated into the 
SLLA as corrections to the intra-level electron- 
electron interaction. 


Note that the perturbatively small corrections to the 
electron-electron interaction in the intermediate regime 
may lead to a qualitative change in the phase diagram. 
Indeed, it is known that small corrections to the pseu¬ 
dopotentials’ values can lead to a transition to a different 
topological phase or to a collapse of the bulk gap (see, 
e.g.. Ref. [2^b 


In this section we discuss factors which determine 
boundaries between the regimes, we also show how to 
take into account the corrections in the third regime. 
First, a brief discussion is presented in subsection III.1[ 
then the technical details are given in subsection |III.2| 


HI.l. Effects important in bilayer graphene 

Now we discuss in details the effects which are impor¬ 
tant in BLG. For the reasons discussed earlier, we con¬ 
tinue suppressing spin quantum numbers of electrons. 

The important effects are as follows: 

• Firstly, BLG in perpendicular electric field is a nar¬ 
row gap semiconductor, thus the effects of vacuum 
polarization are strong [38] . 

• Secondly, Coulomb interaction of electrons can lead 
to mixing of Landau levels, or to emergence of spin 
or/and valley unpolarized states. Coulomb inter¬ 
action can also lead to intervalley hopping of elec¬ 
trons. Even though such processes are suppressed 
compared to intravalley scattering one should still 
estimate their relevance. 

• Thirdly, even when LL mixing is small, virtual 
hopping between the LLs can change (renormalize) 
intra-level electron-electron interaction. 

Next we discuss each of those effects in detail. 



Figure 2: Feynman diagrams showing renormalization 
of the electron-electron interaction due to (a) the vac¬ 
uum polarization processes, and (b) the simplest processes 
involving virtual hopping of one or both of the two interact¬ 
ing electrons from the n-th LL to the n'-ih LL. 


III. 1 a. Vacuum polarization 


The virtual processes shown in Fig. lead to renor¬ 
malization of the electron-electron interaction. This is 
important since the interaction determines the FQHE. 
The Fourier transform of the renormalized (screened) in¬ 
teraction potential can be expressed as^ 


^scrio) 


r(g) 

1 + PV{q)Il{q,aj = 0) 


( 66 ) 


where V{q) = 27Te^/{lqK,) is the Fourier transform of 
the unscreened Coulomb potential, k is the dielectric 
constant, which is felt by the system’s electrons^^, and 
n(g,co’) is the polarization function. Since we are inter¬ 
ested in the effects at the energy scales much less than the 
inter-LL distances we can neglect the retardation effects 
(use only cj = 0). 

We compute the polarization function for the BLG in 
magnetic field within the RPA (random phase approxi¬ 
mation), which can be justified within the 1/A-expansion 
[38] {N = 2 spin projections x 2 valleys = 4). Since 
Il{q^uj = 0) cx: screening is not efficient at large dis¬ 

tances; however, it strongly affects the first few Haldane 
pseudopotentials (corresponding to distances of the order 
of the magnetic length) which have the most significant 
impact on the stability of any FQHE state. Details of 
the calculation are described in sub-subsection [HL2a[ 


III.lb. Landau level mixing and population reversal 

The order of levels in Fig. [^prescribes the natural or¬ 
der of filling of the LLs by electrons in the independent 
electron approximation. However, it can happen that 
for some filling fractions the electron-electron interaction 
leads to a reversal of this natural order in a part of param¬ 
eter space (by external parameters we mean the magnetic 
field, the mini-gap and the dielectric constant). For ex¬ 
ample, the Coulomb energy of the fully filled (2, +1) LL 


^ A similar screening approach has been discussed in GaAs; see, 
for example, na. 

The dielectric constant sensed by BLG is n = +Av 2 )/ 2 , where 

K,i and K ,2 are the dielectric constants of the environment below 
and above the sheet of BLG. 
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is less than the one of the fully filled (2,-1) for U > 0. 
In the region where the interaction is strong compared 
to the gap between the two levels this can lead to the 
fully filled (2, +1) LL having lower total energy than the 
fully filled (2, —1) level. Thus, the former would be filled 
before the latter. 

Whereas for fractional filling such effects are much 
more difficult to analyze, population reversal at the in¬ 
teger filling fraction would be an indicator of a strong 
violation of the SLLA. Thus, we constrain our analysis 
to the region of the parameter space where no popula¬ 
tion reversal occurs at integer filling. More details on the 
population reversal issues are presented in sub-subsection 

inr^ 

When the quasidegenerate levels are from different val¬ 
leys, valley-unpolarized states may be preferred, par¬ 
ticularly for fractional filling. Furthermore, when the 
quasidegenerate levels are from the same valley (as in the 
n = 0 and n = 1 case) level mixing may occur. These are 
interesting effects which are, however, beyond the scope 
of this work. 

In any of the cases (valley unpolarized state or level 
mixing) one has to consider several LLs simultaneously. 
This is hard technically since for the same number of elec¬ 
trons the system’s Hilbert space is significantly larger, 
which complicates use of the numerical diagonalization. 
Therefore, we restrict study to the region of the parame¬ 
ter space where no valley unpolarized states and no level 
mixing occur. Our criteria for smallness of level mixing 
and valley unpolarization are discussed in sub-subsection 
IIII.2bl 

For the valley-polarized states one can still investi¬ 
gate whether the state is spin polarized. Generally, spin- 
unpolarized states are not favored by Coulomb repulsion 
unless the pseudopotential is hollow core (pseudopoten¬ 
tial do not fall off monotonically with relative angular 
momentum). We find that the screened pseudopotential 
does fall off monotonically for > 10 in all the cases con¬ 
sidered in section IV For /^ < 10 the pseudopotential is 
weakly non-monotonic in some regions of the parameter 
space. We, therefore, consider only spin-polarized states 
without further investigation of the restrictions that this 
condition imposes on the parameter space. 


We restrict the region of validity of our consideration 
by requiring the typical interaction energy scale (it can 
be interaction potential value at the magnetic length dis¬ 
tance or, almost equivalently, the zeroth pseudopotential) 
to be smaller than the distances to each of the neighbour¬ 
ing LLs from the same valley. In these regions we take 
into account corrections to the pseudopotentials up to 
second order perturbation theory (Fig. |^). 

The second order perturbation theory also gives rise to 
three-particle interaction [HHni, which can play an im¬ 
portant role in some cases. For example, the three-body 
interaction is crucial to distinguish between the Moore- 
Read Pfafhan and anti-Pfafhan. However, in the present 
work we neglect these terms. 

The details of the calculation of the corrections are 


presented in sub-subsection HI.2c 


III.Id. General plan for numerieal study of a fraetional 
QHE in hilayer graphene 

With the remarks made above, the general plan for 
study of FQHE at a certain filling fraction on a certain 
LL can be formulated as follows: 

1. Calculate the screened interaction potential. 

2 . Determine the region of parameter space in which 
no valley unpolarized states emerge and level mix¬ 
ing doesn’t take place^^. 

3. Calculate pseudopotentials in this region of param¬ 
eter space. Take the corrections due to virtual hop¬ 
ping into account (the modified SLLA). 

4. Use the calculated corrected pseudopotentials for 
numerical diagonalization and compare the exact 
numerically found ground state with the trial one. 

The next subsection contains details of calculation of 
the polarization function, of the corrections to the pseu¬ 
dopotentials, of the criterion for absence of population 
reversal of LLs, and of the criteria for absence of valley 
unpolarization and LL mixing. 


Ill.le. Renormalization of pseudopotentials due to virtual 
hopping 

The SLLA is exact in the limit of infinite energy differ¬ 
ence AE between the LLs. For finite AE the pseudopo¬ 
tentials acquire corrections due to virtual transitions be¬ 
tween the LLs such as, for example, shown in Fig. It- 
Such corrections are theoretically tractable only in the 
perturbative regime (when they are small); however, even 
the presence of small corrections may dramatically affect 
the phase diagram due to the extreme sensitivity of the 
FQHE states to the details of the interaction (see e.g. 

m)- 


III.2. Effects important in bilayer graphene: 
calculation details 

III. 2a. Caleulation of the polarization funetion 

Here we calculate the vacuum polarization function. 
We first do this for the case of integer filling. At the end 
of the sub-subsection we generalize this calculation to the 
case of a fractional filling factor. 


By external parameters we mean the magnetic field, the mini-gap 
parameter and the dielectric constant. 
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The polarization function within the RPA is just a 
density-density correlation function^^ in the free theory 
(this corresponds to the fermionic loop in Fig. [^) 

n(r-r',t-t') = -i(Tp(r,t)p(r',t')), (67) 


h is put to be 1 in this sub-subsection, the T-symbol 
denotes time ordering: 


Tpir,t)p{r\t') = { J f 

Density is meant to be nbrmW braced? electrbn/hole creation 
operators should be to the left of the annihilation operators. 


The density-density correlator is translation-invariant 
since the system is uniform, so 

{Tp{r,t)p{r\t')) = {Tp{r-r\t-t')p{0,0)). (69) 


Let’s denote the set of quantum numbers (n, m, s) by 
k, and write k < kp the state is occupied, and k > kp 
otherwise. The polarization function can be expressed 
with the help of the wave functions (47) in the following 
way: 


n(r, t) = — 2i 


k<kF,k'>kF,C=^' 


x 

k^kF fk' <CkF 


, (70) 


0(x) = 


1, X > 0 
0, X < 0 


(71) 


The factor of 2 in front of the square brackets is due to spin. 

The Fourier transform of the polarization function is then defined as 




(72) 


So the polarization function at the zero frequency n(q, cj = 0), which we need to find the interaction potential, can 
be expressed as 


n(q,w = 0 ) = 2 ^ TP ^ TP *‘*''/q 5 ffc(x)^^'fc/(a:;) 5 ffc/( 0 )^^'fe( 0 ) + c.c.) . 

7^7 7 7/7 7- r-/ -^k' J 

k>kF ,k' KkF 


After a short calculation we find that 

2 


2 _^ 2 1 

II(q, C(J ~ ~ J2 ^ ^ ^ {In,s,n' ,s' T C.C.) = 1‘^fyjj ^himless(^)5 




-t, 


(73) 


(74) 


(_!)(« "Oy drrJo{qr/l)x V’n,o(w')V’n',o(w) + 

WBiU^_2,2{w)i’n'-2,2{w) + (t^)) + 

(^AiUn,-2Mi’n',-2{w) + BiUn-2,oMi’n'-2,oH + 

)TlPn-l,-l{w)TP^,-i,-i(w)) + 

(C^®C|f' + lpn-liw)^n' -i{w) + ^n- 2 ,l{w)^n'- 2 ,l{w) 

+ (C|*C^f + )7/)„_i,o(u;)V’n'-l,o(w 


(75) 
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where j denotes a set (n,^,5), overbar denotes complex conjugation^^ The integrals in the expression above can be 
found analytically: 


f 


dr rJo{qr/l)'ipn,m{w)'ip„ 




hi 




i(w^ — o 

Zir 




(76) 


i2 

-^0,^2(^)In>*2 ~ 

k=0 


k\{k + ii - i2)K^2 - k)\ 





ii-i2 


T (n-n) 
^i2 


{\Am. 


(77) 


■ « 1,«2 


(^)|ii<i2 — -^*2,0 (^) 


(78) 


where x = Qx ^ iQy^ and lA are generalized Laguerre polynomials. The derivation is very similar to the derivation 
of formula (32) presented in Appendix [A| 


The Fourier transform of the interaction potential (31), 


(66) can be expressed then as 


^scrio) — 


27re^ 1 

nql 1 _L e2 Hd imless (q) 


(79) 


In the case of a fractional filling factor we have to take 
into account not only fully filled or entirely empty LLs 
but the partially filled ones as well. We do this with the 
help of the following approximation: for a partially filled 
level we add terms which correspond to the level as an 
empty one and as a filled one with coefficients (1 — {u}) 
and {z/} respectively. For example, if some level is half- 
filled then all the terms in the polarization function which 
correspond to the hopping to this level appear with the 
coefficient 1/2 = 1 — 1/2, and the terms which correspond 
to hopping from this level also appear with the coefficient 
1/2. Thus, we do not take into account correlations in a 
partially filled LL. 

For this work the polarization function was calculated 
approximately: we calculated only terms with n,n' < 
^cutoff = 4. We checked that the pseudopotentials in the 
region we are interested in differ from the pseudopoten¬ 
tials calculated with Ucutoff = 3 by less than 2%. 


III.2b. Population reversal of Landau levels and level 
mixing 

It was discussed in sub-subsection lTIl. Ibi that when the 
typical energy scale of the electron-electron interaction 
becomes larger than the difference of kinetic energies of 
two LLs from the same valley it is natural to expect the 
SLLA to break down. The numerical study in such cases 
is significantly hampered. Moreover, it is hardly probable 
to find a state similar to a single-level state in this regime 
of strong level mixing. Therefore, we would like to work 

Note that for some values of n, n' this general expression has 
the non-defined wave functions like, e.g., —2 o('^) ^^r n' = 0 

OT n' = 1. These terms, however, do not contribute, which is 
ensured by the coefficients 5^, , Cn^ etc. which take zero values 
in those cases. 


only in the regime where the mixing of LLs is small. For 
this we demand the typical energy scale of the electron- 
electron interaction to be smaller than the kinetic energy 
distance to the closest LL from the same valley. The 
remnants of the level mixing can be incorporated then 
into the corrections to the SLLA which are discussed in 
the next sub-subsect ion. 

One can use different quantities to define the typical 
energy scale of the electron-electron interaction. For ex¬ 
ample, one can use the interaction potential value at the 
magnetic length distance V{1) or the zeroth pseudopoten¬ 
tial at the LL one is interested in They typically 

differ by a factor of order of unity, which is not too im¬ 
portant. We choose to use value of the zeroth pseudopo¬ 
tential as the typical interaction energy scale. Therefore, 
we restrict the region of consideration to those values of 
external parameters 17, 5, n for which < AF^, 

where /S.E is the kinetic energy difference between the 
LL under consideration and the closest other LL to it. 

There is a subtlety regarding this restriction. One can 
use the zeroth pseudopotential for the screened or for the 
bare Coulomb interaction potential. Using the Coulomb 
pseudopotential seems natural as it is the fundamental 
perturbation theory controlling parameter. However, for 
the weakened screened potential the Landau level mixing 
is smaller, and restricting the applicability region by the 
Coulomb interaction scale one excludes regions where our 
approach should still give reliable results. On the other 
hand, if the screening is so strong that the Coulomb en¬ 
ergy scale significantly differs from the screened one, then 
the RPA approach we use for calculation of the screened 
potential may be not good enough, bringing in an error 
in the interaction potential. In section [TV] we restrict our 
region of consideration by the screened energy scale, how¬ 
ever, for comparison we also show the region’s boundary 
calculated for the unscreened Coulomb potential. 
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Apart from mixing of LLs from one valley, a similar 
process can take place for neighbouring levels from dif¬ 
ferent valleys like (2,±1) (see Fig. [^. This is due to 
Coulomb interaction on the lattice scale that can make 
electrons jump from one valley to another. This inter¬ 
action is considered in more detail in [39]. What is im¬ 
portant for us is that the typical energy scale of this 
interaction is 


^^Coul(^) 



27re^ 3e^a 


(80) 


where a ^ 0.25nm is the graphene lattice constant. We 
would like this interaction not to play a significant role. 
Therefore we restrict the parameters region by demand¬ 
ing that its energy scale is smaller than the distance be¬ 
tween the LL under consideration and its closest neighbor 
from the different valley. 

Even if the valley index is a good quantum number 
(i.e. electron jumping between valleys is suppressed), 
the ground state might not be valley polarized. The rea¬ 
son for this can be the competition between the kinetic 
energy favoring a LL from one valley and the exchange 
energy favoring a LL from the opposite valley. If this hap¬ 
pens, the two LLs, even not mixing, influence each other 
through the density-density interaction (since the elec¬ 
trons still repel each other). In such case the two levels 
from different valleys should be considered together just 
like in the case of level mixing. So by the same reasoning 
as in the case of mixing, we do not want to consider the 
system in the regime of two levels from different valleys 
partially filled. 

Thus, we need to find the region of parameter space 
where such simultaneous filling does not occur — in order 
to use the SLLA there. Since the case of a partially 
filled level is hard to analyze, as a criterion we choose 
to demand that for integer fillings there should be no 
change of the filling order. I.e., the full energy (kinetic 
plus interaction) per electron of fully filled levels should 
put them in the same order as their kinetic energy. For 
example, if the kinetic energy of the (2,-1) LL is less 
than the one of the (2, +1) LL, then the full energy per 
electron in the fully filled (2,-1) level should also be less 
than the full energy per electron in the fully filled (2, +1) 
level. 

For the fully filled level it is easy to calculate its in¬ 
teraction energy since there is only one state possible — 
the Slater determinant of all the level’s orbitals. The 


Unlike the case of level mixing in one valley, due to quasi¬ 
momentum conservation, mixing of the LLs from different valleys 
can happen only if both of them are filled with electrons at least 
partially. Naively, one would think that because of this argument 
the level (2,-hi) is not dangerous when we consider the (2, —1) 
LL. However, the screening processes happen because of hopping 
of the electrons to higher LLs. Therefore, the (2, -hi) LL is ’’vir¬ 
tually” filled, to some extent. Thus, to be on the safe side, we 
still apply this restriction when we consider the (2, —1) LL. 


interaction energy can be expressed through the pseu¬ 
dopotentials: 


L/inter., N electrons — Vp2, (81) 


E; 


inter, pp 


lim 

N^oo 


L/inter., ^electrons 

N 


2 F 

me2N-l 

(82) 


where p 2 is the two-particle density matrix, pp in the 
subscript stands for ’’per particle”. 

This energy can be separated into the Hartree (density- 
density interaction) and the Fock (exchange) contribu¬ 
tions: 


-^Hartree 


PP 


oo 

E 

m=0 




1 



dr rV{r), 


(83) 


^Fock pp = ^(-1)’”+T„ 
m=0 


-T 

Jo 


dr rV{r)g{r). 


(84) 

The function g{r) is related to two-particle and one- 
particle density matrices p 2 and pi\ 

pi{x’\x) = {x’\pi\x), (85) 

p 2 {x'y'\xy) = (a;'| 0 {y'\ p 2 h) 0 jy), (86) 


P 2 {x'y'\xy) = 

{pi{x'\x)pi{y'\y) - pi{y'\x)pi{x'\y)) , (87) 


g{r) = Ar^(p 2 (?’, 0 |r, 0 ) - /?i(0|0)^) = 

-7VVi(0|r)pi(r|0). (88) 

Note, that the Hartree energy can be expressed as an 
integral of the interaction potential, with the form of the 
integral independent of a Landau level. This is a conse¬ 
quence of the fact that the fully filled LL has constant 
density. We note that the integral for the Hartree en¬ 
ergy is divergent at the upper limit for the Coulomb-like 
interaction potentials. However, only differences of the 
energies have physical meaning, thus we can calculate 
this integral with a certain regularization if we keep the 
regularization always the same. 

The Fock contribution is, on the contrary, convergent, 
but it depends on the LL through the function ^(r), 
which characterizes short-range correlations. 

If the interaction potential V (r) is the same for two dif¬ 
ferent Landau levels, their Hartree energies are identical. 
However, as the screening can be different for different 
LLs, their Hartree energies (83) are different. One may 


think that this energy difference can contribute to pop¬ 
ulation order reversal. However, one has to take into 
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account the background positive charge (since the sys¬ 
tem is electrostatically neutral). The total electrostatic 
energy then, as we show in Appendix does not differ 
for differently screened potentials. 

Therefore, the interaction energy difference comes from 
the Fock term only. For the non-relativistic levels and the 
Coulomb potential, the magnitude of the Fock energy de¬ 
creases with the level number n, while having a negative 
sign. Thus, both the Coulomb interaction and the kinetic 
energy favour the natural LL population order. This is 
not the case in bilayer graphene. 

Consider, for example, two levels (2,-1) and (2,-1-1) 
in BLG. The latter level has a greater kinetic energy. 
The wave functions in the (2, +1) level are close to the 
wave functions of the non-relativistic n = 0 LL, while 
the wave functions in (2,-1) are close to the ones in 
n = 2. Thus, the Fock energy prefers the (2, +1) LL, 
while the kinetic energy prefers the (2,-1) LL. Therefore, 
if the interaction is strong enough population reversal 
may happen. 

While for the bare Coulomb interaction population re¬ 
versal would happen in some regions of the parameter 
space, for the screened potentials we do not find such 
an effect for the (2,-1) and (2,+1) LLs. This makes 
emergence of valley-unpolarized states improbable. 

In summary, in this sub-subsection we considered 
the effects of mixing of LLs from one valley, mixing 
of LLs from different valleys, and emergence of valley- 
unpolarized states. We find that the last one does not 
occur in realistic conditions. The other two effects re¬ 
strict the applicability of the SLLA. 


troduce the following basis in the Hilbert space: 


{ |m,M,jf,jf)), if jf = j', 

^ , II J ^ J , 

(89) 

where m G 2Z+ -h 1, M G Z+, and 


\m,MJj')) = 


V 2 ^+^m!M! 

{b\ - birib\+ (90) 


Note that the basis vectors ( [8^ are antisymmetric with 
respect to electron permutations, while the auxiliary vec¬ 
tors ( [^ are not. 

The leading correction to the eigenstates’ energies, 
which is due to the virtual hopping to the higher LLs 
from the same valley, is given by the second order per¬ 
turbation theory: 


=V^- 

m m 

E 

{ji,j2<ji)9^{j,j),m' ,M' 


\{m,M,j,j\V\m',M',ji,j 2 )\" 


^kin _|_ ^kin _ 2^^kir] 


(91) 


Since the interaction potential F is a function of 
one can show that the only non-zero matrix elements of 
all the {m,M,j,j\V\m',M',ji,j 2 ) are 


{m,M,j,j\V\m+ (ni - n) + (n 2 - n), M, ji, j 2 ) = 

(92) 




IIL2c. Calculation of corrections to the SLLA 
pseudopotentials due to virtual hopping 


Suppose we are in the region where spin-/valley- 
unpolarized states do not emerge, and the LL mixing 
of the level under consideration with the levels from the 
same valley is small. Then the small mixing can be taken 
into account with the help of corrections to the electron- 
electron interaction within the SLLA. Those are the small 
corrections to the pseudopotentials. It is known that 
small corrections of the order of 5 — 10% to the pseu¬ 
dopotentials’ values can lead to a significant change of 
the overlap with a trial state (see, e.g.. Ref [23]). In this 
sub-subsection we present the formulae we use to com¬ 
pute such corrections. 

Consider the two-particle problem. In the subsection 
im it was shown that the eigenstates of the two-particle 
problem within the SLLA are |m, M), with their ener¬ 
gies being j = Let us denote |m,M) as 

to emphasize that both of the electrons are 
in the LL j. Now we add to our consideration the closest 
unfilled LLs from the same valley (for small deviations 
from the naive SLLA those are the levels above); we in- 


yj,3,01,32 _ 

^ m 

{{m,M,j,j\V\m+ -n) + (n 2 - n),M, ji, ja)), 

(93) 

where 5j^j^ is equal to 1 when = j 2 and 0 otherwise. 
We have used the fact that 

y3,3,31,32 _ (^_Y^^{ni-n) + {n2-n)yj,j,j2Si ^ 

The auxiliary matrix elements can be ex¬ 

pressed through non-relativistic matrix elements, simi¬ 
larly to how the pseudopotentials in BLG are expressed 
via the non-relativistic pseudopotentials: 


|m, M, n, n')) 


{b\ - iyA + by -n')2, (95) 


yn,n ,ni,n 2 _ 

* m 

{{m,M,n,n'\V\'i 
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yrn,n,ni,n2 _ /_i\(ni—n) + (n 2 —n)-^n,n,n 2 ,ni 


(97) 


yodJi ,32 

^ m 



T/n-2,n,ni-2,n2 td^s jd^si a^S2 i 
Vm ^n2 


(98) 


The non-relativistic matrix elements can be expressed 
through the Fourier transform of the pseudopotential in 
a form quite similar to the expression for the pseudopo¬ 
tentials (32): 


V 

V r 


n,n ,ni,n2 _ 


poo 

I —n)-|-(n 2 —n') (^V^) X 

Jo 

F„,„,(a:)F„,,„,(-x)e-«'^, (99) 

Ztt 
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|a;P\^ f ix V 
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@ ( 100 ) 


-^ 21,22 (^) |il <22 — -^22 


l(^). 


( 101 ) 


X = Qx ^ '^Qy^ and are generalized Laguerre poly¬ 
nomials. The derivation is very similar to the deriva¬ 
tion of formula (32) presented in Appendix H Since 
Fr, n{x) = Ln{\x\y2)^ for n = ni, n' = n 2 the formula 

S is restored, as it should be because by definition 

15212) _ -^ni,712,211,712 

y jYi 

There is a subtlety regarding what interaction poten¬ 
tial V should be used in the formulas above. It is tempt¬ 
ing to use the screened interaction potential (66) to calcu¬ 


late virtual hopping corrections. However, as we show in 
Appendix this would exceed the accuracy of the per¬ 
turbative calculation. Therefore, we use the unscreened 
Coulomb interaction potential in calculating the virtual 
hopping corrections to the SLLA pseudopotentials. 


III.3. Summary of the section 

In this section we discuss the effects which restrict the 
applicability of the SLLA in BLG. We also discuss the 
conditions under which it is sufficient to introduce cor¬ 
rections to the SLLA to restore the theory applicability. 
We calculate these corrections in the second subsection. 


IV. POSSIBILITY TO OBSERVE THE 
MOORE-READ STATE IN BILAYER GRAPHENE 

In order to investigate the role of the effects discussed 
above on the stability of FQHE states we focus on the 


Moore-Read Pfafhan. Our choice is motivated by the 
following considerations. Firstly, this state is particularly 
sensitive to the details of the interaction so it is a good 
illustration for our analysis. Secondly, the stability of 
this state in BLG was investigated in Refs. Eiiini in the 
SLLA approximation but without these effects taken into 
account, so we can compare the phase diagrams. Thirdly, 
the Moore-Read state itself is an important state because 
it is an example of the non-abelian topological fluid. 

The tunable parameters are the magnetic field H, the 
electric field which determines the mini-gap parameter 
U and the effective dielectric constant which controls 
the deviation from the naive SLLA (which is exact for 
K ^ 00 ). We can also choose the half-filled LL num¬ 
ber. Here we will concentrate only on the two levels: 
(1,-1) and (2,-1). The (1,-1) level wave function is 
constructed from the nonrelativistic n = 0 and n = 1 LL 
wave functions, the (2,-1) level wave function is con¬ 
structed from the nonrelativistic n = 0,1,2 LL wave 
functions. In both cases for the bare Coulomb interaction 
one can tune the pseudopotentials close to their values 
at the nonrelativistic n = 1 LL, where the 5/2 state is 
observed in GaAs^^. 


The tuning mechanisms are, however, different for the 
two levels. Amplitudes of the wave function (47) in the 
(1,-1) LL show little dependence on U so the main con¬ 
trol parameter is . In contrast, the amplitudes of the 
wave function in the (2,-1) LL mainly depend on one 
parameter which is the U/hWc ratio, so both B and U 
can be used for tuning. 

The main factors determining deviation from the naive 
SLLA for the two levels are the polarization and virtual 
hopping to the nearby levels. For the (1,-1) LL this is 
hopping to the (0, —1) and the (2, —1) LLs, while for the 
(2,-1) LL the important hopping is to the (3,-1) LL. 
In addition to this, for the (2,-1) LL, it is important 
to consider effects of mixing with the (0, —1) and (2, +1) 
LLs. The latter are important factors restricting the re¬ 
gion of applicability of perturbative analysis, however, 
when suppressed they do not lead to a renormalization 
of the intra LL interaction. 


Figures and show the regions of the applica¬ 
bility of perturbative analysis for different values of n 
for the (I, —I) LL.^^ For Fig. the typical interaction 


The dielectric constant sensed by BLG is k = (ki +/C 2 )/ 2 , where 
Ki and K ,2 are the dielectric constants of the environment below 
and above the sheet of BLG. 

Though the (2, +1) level wave function is also constructed from 
the nonrelativistic n = 0,1,2 LL wave functions, the numerics 
shows that the high overlap with the Moore-Read state is not 
achieved here for the bare Goulomb interaction in contrast to 
the (1, —1) and the ( 2 , — 1 ) LLs. 

In the low-energy two-band model HU (which corresponds to the 
71 —)■ 00 limit) such tuning is impossible because the amplitude 
is identically equal to 1 with other amplitudes being zero. 
Due to some errors in calculation of the region of applicability, 
in the original result-reporting paper [?] the form of the region 
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Figure 3: The region of the applicability of perturba¬ 
tive analysis for fixed values of = 5 (yellow region), 
10 (green and yellow regions) and 15 (blue, green, 
and yellow regions) for the (1,-1) LL. The size of the 
region increases with increasing k. For (a) the typical inter¬ 
action energy scale is taken to be the zeroth pseudopotential 
of the screened interaction potential, for (b) — of the bare 
Coulomb potential. The thick black line shows where the 
maximum overlap with the Moore-Read Pfaffian state for the 
bare Coulomb interaction is achieved. (Color online). 


energy scale, which determines the significance of level 
mixing, is estimated with the help of the screened poten¬ 
tial (’’type S estimate”), while for Fig. — with the 
help of the bare Coulomb potential (’’type C estimate”). 
The regions are bounded from above by the condition of 
small hopping to the (2,-1) LL, the lower bound is due 
to the condition of small hopping to the (0, —1) LL. At 
small enough magnetic fields at least one of the condi¬ 
tions is violated at all values of U. The thick black line 
shows where the maximum overlap^^ with the Moore- 
Read Pfaffian for the bare Coulomb interaction (which 
is about 0.95, compare with non-relativistic n = 1 level 
overlap of 0.7) is achieved. One can see that for small 
dielectric constants this line lies outside the region of va¬ 
lidity of perturbative analysis, however for large enough 
K they intersect near U = bO meV. 

Figures and[^ show the regions of the applicability 
of perturbative analysis for different values of k for the 
(2, —1) LL. For Fig.[^ the type S estimate is used, while 
for Fig. S' the type C estimate is used. The regions are 
bounded from above by the condition of small mixing 
with the (0,-1) LL, the lower bound is due to the condi¬ 
tion of small mixing with the (2, +1) LL. The thick black 
line shows where the maximum overlap with the Moore- 
Read Pfaffian for the bare Coulomb interaction (which 
is about 0.92, compare with non-relativistic n = 1 level 
overlap of 0.7) is achieved. One can see that for small 
dielectric constants this line lies outside the region of va¬ 
lidity of perturbative analysis, however for large enough 
hz they intersect near U = 30 meV. 

Figures and|^ show the dependence of the overlap 
of the exact ground state of the system with the Pfaf¬ 
fian on the magnetic field and the dielectric constant at 
U = bO meV for the (1,-1) LL and at = 30 meV 


Figure 4: The region of the applicability of perturba¬ 
tive analysis for fixed values of = 2.5 (brown region), 
5 (yellow and brown regions) and 10 (green, yellow, 
and brown regions) for the (2,-1) LL. The size of the 
region increases with increasing k. For (a) the typical inter¬ 
action energy scale is taken to be the zeroth pseudopotential 
of the screened interaction potential, for (b) — of the bare 
Coulomb potential. At At = 2.5 in (5) the condition of small¬ 
ness of level mixing is not satisfied anywhere within the range 
of external parameters shown. The thick black line shows 
where the maximum overlap with the Moore-Read Pfaffian 
state for the bare Coulomb interaction is achieved. (Color 
online). 
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Figure 5: Color plot of the overlap of the ground state 
with the Moore-Read Pfaffian for 12 particles as a 
function of the magnetic field B and the dielectric 
constant K. (a) - for the (1,-1) LL at f/ = 50 meV, (b) - 
for the (2,-1) LL at = 30 meV. Contours show the lines 
of constant overlap. The region where perturbative analysis 
is not applicable according to the type C estimate is hatched. 
Data is not shown beyond the region where perturbative anal¬ 
ysis is applicable according to the type S estimate. (Color 
online). 


for the (2,-1) LL respectively.^^ We do not show the 
data in the region where the perturbative analysis is not 
applicable according to type S estimate. The region of 
inapplicability of our theory according to the type C esti¬ 
mate is hatched. As one can see, for the (1,-1) level, the 
highest overlap achieved at dielectric constants as high as 
hz = 40 is about 0.7 despite having overlaps up to 0.95 for 
the bare Coulomb interaction (which corresponds to the 
K ^ oo limit). For the (2,-1) level a high overlap up to 
0.92 is achieved at k ^20. 

The authors of [6] found that in the (1,-1) LL, high 


is not entirely correct. This applies to both the (1, —1) LL and 
the (2, —1) LL case, which is considered next. 

By overlap we mean scalar product’s absolute value squared. 


Due to some mistakes in calculation of the screened potential, in 
the original result-reporting paper [7] the overlaps for the (1,-1) 
LL behave somewhat differently. Here we corrected the mistakes. 
As to the (2,-1) LL, here we consider another value of U than 
in the paper [7], see footnote |18| 
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overlap is achieved in the region near B = 10 T. We 
find that the region of high overlap is situated there for 
K ^ 60. However, for realistic values of k the effects 
of level mixing become significant which makes obser¬ 
vation of the Moore-Read state unlikely, especially for 
K < 10.^^ The (2,-1) LL was also considered in [6], 
where it was concluded that the maximal overlap with 
the Moore-Read Pfafhan on this level is less than 0.6. 
Our results do not support this conclusion (even for the 
bare Coulomb interaction). 

Thus, our results show that the (2,-1) LL is a bet¬ 
ter candidate to observe the Moore-Read Pfafhan state 
in BLG than the (1, —1) LL. However, even in the case 
of the (2,-1) LL to tune into the regime of high over¬ 
lap with the state one needs a dielectric constant about 
hv = 20. This is much higher than the usual k ^ 2.5 
for graphene on Si02 or hexagonal boron nitride (hBN) 
substrate. Even on Hf02 substrate [42] k is around 12.5. 
A possible way to achieve higher dielectric constants is 
to use substrates on the both sides of the BLG sheet. In 
that case the effective dielectric constant k is equal to 
the substrate dielectric constant. Therefore, using Hf02 
one can get n = 25. 

The region of high overlap intersects this value of di¬ 
electric constant for the (2,-1) LL. For example, at 
U = 30 meV, B = S T the overlap of the exact ground 
state and the Moore-Read state is about 0.92, with the 
gap to the first excited state being about 1.6 K. With 
increasing magnetic field the gap monotonically increases 
to the values of around 5KatH = 15T. At the same 
time the overlap decreases to around 0.7 which is still 
fairly large. Similarly, with decreasing the dielectric con¬ 
stant down to /^ = 10 (for which one still needs Hf 02 
substrate but only on one side of the BLG sheet) the gap 
monotonically increases to the values of around 3 K, with 
the overlap decreasing to 0.4 which is not too small. 
This result, obtained for a finite number of particles, sug¬ 
gests that the system may still be in the same topological 
phase at higher magnetic fields and lower dielectric con¬ 
stants. 

However, the experimental observation of the Moore- 
Read Pfafhan in BLG in this way in the near future is 
unlikely. The main problem for observation of the frac- 


It is interesting to note that recently a = 1/2 FQHE state has 
been observed in suspended bilayer graphene [9] for U ^ 0. The 
parameters of the experiment he far outside the region of appli¬ 
cability of our methodology. However, the study m, which does 
take into account the strong mixing between the quasidegener¬ 
ate (1,-1) and (0, —1) LLs beyond perturbation theory, claims 
that the Moore-Read Pfafhan is a likely candidate to explain the 
observed u = 1/2 FQHE. 

Such conhgurations with hBN as a substrate have recently 
started being explored experimentally from the perspective of 
FQHE: see Ref. [To] . 

In our simulations the Hilbert space is about 16000-dimensional. 
Therefore, two random vectors would have the average overlap 
about 1/16000. 


tional QHE in graphene (and BLG as well) is the too high 
disorder in the samples caused by the disorder in the sub¬ 
strate (evidently, it is impossible to observe FQHE when 
the typical height of the disorder potential is greater than 
the gap to the first excited state). For example, FQHE 
with the filling factor = 1/3 has been observed in sus¬ 
pended single-layer graphene |3| and in graphene on hBN 
substrate [5| but not on a substrate with a different lattice 
structure (which is the case for Hf 02 needed to achieve 
the high dielectric constant). 

The reader can find some additional data on the be¬ 
haviour of gaps and overlaps in Appendix [D| 

One important thing to note in the context of modern 
experiments is that the possibility to tune the mini-gap 
parameter U is achieved through placing two metallic 
gates on the both sides of the BLG sheet, very close to it 
m- Therefore, one can expect that the screening of the 
Goulomb potential will be due to not just the internal 
BLG dynamics but also due to the metallic gates. We 
have not taken this effect into account in our analysis. 


Conclusions 

We analyze the influence of inter-Landau level tran¬ 
sitions (vacuum polarization, virtual hopping) on the 
phase diagram of the FQHE states. We find that the 
SLLA can only be used under quite stringent conditions, 
and corrections to the SLLA should be taken into ac¬ 
count. A method for taking the corrections into account 
by means of perturbation theory is developed. However, 
the region of applicability of the SLLA with our correc¬ 
tions is also quite restricted. 

With the help of the developed method we study the 
possibility to observe the Moore-Read state in bilayer 
graphene. We find that the mentioned effects indeed lead 
to a substantial modification of the phase diagram. How¬ 
ever, the external parameters needed to tune into the 
regime favouring the Moore-Read state are, in principle, 
achievable, for BLG surrounded with Hf02. 

We also developed a set of programs in Wolfram Math- 
ematica which implement our methodology to study the 
case of Moore-Read Pfafhan state. To obtain the pro¬ 
grams please contact K.S. 
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Appendix A: Derivation of the expression for 
pseudopotentials in terms of the Fourier transform 
of the interaction potential 

In this appendix we derive an integral representation 
for the pseudopotential 

= (ni,n2,m,M|l/|ni,n2,m,M) = 

(ni,n 2 ,m, 0|y|ni,n2,m,0), (Al) 
where the two-electron states 

Ini, 712, m,M) = , ^ X 

ib\ - blr{b\ + (A2) 


have been defined in section |L2| Eq. ( |29| ). 

The direct Fourier transform for the interaction po¬ 
tential is defined in (31). The inverse Fourier transform 
looks like 


Vir) = I 


(Pq 


V{q\ 


,iqr/i 


(A3) 


We can rewrite the definition as follows: 

= (ni,n2,m,M|l/|ni,n2,m,M) = 

/ 72 

(g)(ni, n 2 , m, ^ 2 , m, M). (A4) 

Introducing x = Qx -\- we can write qr// = {xw -h 
xw){2. re, w can be expressed through the d, d’*', 6 , 6 ’!' op¬ 
erators: 


w = \/ 2 (d + 6 ^), w = V2{d) -h h). 


(AS) 


where 




I’^i) = 

di|(l) = 


1 


Thus, the matrix element can factorized into a product 
of three different matrix elements: 

(ni,n2,TO,M|e*‘i“'/'|ni>^2,TO,M) = Mn^Mn^Mm, 

(A 6 ) 

(A7) 
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V^ 2 ! 
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( 61 - 62 ) 111 ) = 


\/2^m\ 

0 . 
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We note that the expressions (A7), (AlO), (A13) are 


similar and can be calculated using the same technique. 
We only give the details of the calculation of Aim- By 
virtue of the Baker-Campbell-Hausdorff identity. 




It follows that 

Mn, = = 


(A16) 


(ni|e7i*“ie7f*“i|ni) = 

00 

Q-\x\ (A17) 


n=0 


where we used the resolution of identity 1 = l^)(^l' 

After a calculation of the expression 
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and its complex conjugate, we get 
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where is the Laguerre polynomial. 

Calculation of Mn 2 ^ Alrn is done the same way. Gath¬ 
ering all the three matrix elements together we finally get 
formula ([3^: 


y{ni,n2) 


poo 

/ r(</)L„(</2)L„^(g2/2)L„,(gV2) 

Jo 


. (A20) 


Appendix B: Full Hartree energy of a Landau level 

The electrostatic (Hartree) energy per electron of a 
fully filled Landau level can be expressed as 

i^Hartree PP = ^ J Vpi^)pivW{\^ “ YI), (B1) 
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where N is the number of electrons in the LL, x and y 
are the coordinate vectors of points in the plane, V{r) 
is the electron-electron interaction potential, and p(x) is 
the electron density, integration is done over the whole 
plane^^. Using the fact that for a fully filled LL p(x) = 
l/(27r/^) = p and changing integration variables to R = 
(x + y)/2 and r = x — y, we get 


such difference through the interaction potential Fourier 
transforms V{q) defined in Eq. (31): 




( 1 ) 

Hartree pp 


-E, 


( 2 ) 

Hartree pp 

— lim 
47r g->o V 


^\q)- 


(B5) 


-S'Hartree pp 


2N 


/ d?RdPrV{r)^ (B2) Recalling (66) and taking into account that n(g) 

U{q,uj = 0) r\j const X + O(g^) for g ^ 0, one gets 


where r = |r|. f d?R is equal to the sample area S. 
p X S = N, therefore. 


-^Hartree pp — ^ j ^ tV — 

If 1 

^27r J drV{r) = ^ dr rU(r), (B3) 


in agreement with Eq. (83). 

This derivation is va! 
sufficiently fast as r 


id for the potentials that decay 
oc, and the Coulomb interaction 
does not satisfy this condition. However, only differences 
of the energies have physical meaning. If we consider 
two energies for differently screened interaction poten¬ 
tials (66): 


Hartree pp Hartree pp 

^ dr r (r) - 1 /( 2 ) ^ (b4) 

the integral of the difference is convergent since the 
Coulomb part cancels out and what remains decays suffi¬ 
ciently fast at large distances. Therefore, we can express 
Over the whole sample of area S' —)■ oo with N electrons in the 
LL so that S/N = 27r/^. I is the magnetic length. 


Hartree pp Hartree pp 



lim 

q^O 


/n(2)(g) -n(i)(g) 
V ^2 


(B6) 


Based on this consideration, naively one might think 
that Hartree energies contribute to the population order 
reversal effect and that the contribution is characterized 
by the quantity lim^^o n(g)/g^. However, this consider¬ 
ation does not take into account the electrostatic energy 
of interaction with the compensating positive charge (as 
the system is electrically neutral as a whole). Although, 
the screening processes happen at the BLG plane, they 
influence the interaction potential of charges outside of 
the plane. 

Consider a simple model: the compensating positive 
charge is evenly distributed in the plane at a distance 
d from the BLG plane; for simplicity, all the system is 
in the environment with the dielectric constant k, = 1. 
Suppose the screening processes happen only in the BLG 
sheet. Then one can calculate the potential of interaction 
of the charges situated in different points of space. We 
are going to do that now. 


We need several definitions: Coulomb electrostatic potential p (^r = -\- y‘^ -\- = Ifr^ the spatial polarization 

function with screening happening only in the z = 0 plane n(r) = Il{x,y^z) = Il{x,y)S{z)^^, and their Eourier 
transforms w.r.t. plane coordinates x^y defined similarly to Eq. (31) and to Eq. ( [7^ respectively: ip{q = |q|, 2 :) = 
27rexp (—g| 2 :|//)/(g/) and Il{q,z) 


U{q)6{z). 

Then the interaction potential of two charges ei and 62 placed in r and r' is 


Uscr(r,r') = eie2Pscr{\^-y\^z,z') = 

eie2p{r,r')- 

- 61626 ^ J d^xd^yp{\T - x|)n(x - y)p{\y - r'|)+ 

+ 61626 ^ J d^xd^yd^zd^wp{\r — x|)n(x — y)p{\y — z|)n(z — w)(/:?(|w — r'|) — ... (B7) 

Here e is the elementary charge — the absolute value of the charge of the electrons participating in the screening 
processes. 


S(z) is the Dirac delta-function. 
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One can show then that the in-plane Fourier transform of the screened interaction potential of the two charges has 
the form 




Vsciiq, Z, z', 61,62) = eie 2 <^scr( 9 , z') = 

e,e20 (,., - + ' 

'2»exp(-:^^) 2irexp(-!iP) eVn{p) 2xexp 


6162 


ql 


ql 


l + e 2 / 2 ^n(^) 


ql 


(B8) 


For ei = 62 = e and z = z' = 0 one restores the expression (31), as it should be. 

The electrostatic energy of interaction of all the electrons in the fully filled LL and the charge-compensating 
background is equal to 


-Bnartree = 6^ j (fx (fyp(x)p{y) ^i(^scr(|x - y\, Z = z' = 0) + 


+ 2 ‘<^scr(|x-y |,2 = z' = d)- (^scr(|x-y|,2: = d,z' = 0)J. (B9) 

Performing the same operations as in the beginning of the current appendix for the Hartree energy of electrons in the 
LL only, we get: 


e^N , f 1 ^ 1 ^ ^ ^ 

-S'Hartree ~ ( q'^ scr(^: Z — Z = 0) + X‘^scr(^: Z — Z = cT) ^scviSli Z> = Z =0) 

ZTI \ Z Z 


2^ (“^)) 47r2e2n(g) 


lim 
27r <7^0 


2q^ (l + 6^P^n{q)'^ 


qd 


1 + exp — 2 ^ —2 exp —- 


qd 


I 


e^N 27rd _ e^Nd 
~ /2 ’ 


(BIO) 


where we have used the fact that n(g) = li^q^uj = 0 ) ^ const x -h 0{q^) for g —)► 0 . 


One can see that the answer does not depend on the 
polarization function 11 (g) and coincides with the energy 
of a capacitor made of two parallel planes of area S at 
a distance d with charge Q = eN on its plates. Indeed, 
such a capacitor has the capacitance C = S/{A'Kd) and 
energy F^cap = I^C). Since S = 2ttPN, one easily 

finds that F/jjartree — -^cap* 


So, in this simple model, the full Hartree energy of 
a filled Landau level and the compensating charge layer 
does not depend on screening. Of course, this model does 
not take into account many features of real experimen¬ 
tal systems. However, it clearly illustrates that one can 
hardly expect to have different Hartree energies for differ¬ 
ently screened potentials. Therefore, the Hartree energy 
does not contribute to the population reversal effect. 


Appendix C: Peculiarities of calculations of virtual 
hopping corrections 

In sub-subsection [hi. 2c| we presented the derivation of 
the corrections to the pseudopotentials which correspond 
to the diagram in Fig. It- The derivation presented there 
is based on a consideration of two interacting particles, 
which interact through a potential, within the method¬ 
ology of first quantized quantum mechanics. In fact, the 
full expression for this correction comes from the con¬ 
sideration within a quantum field theory methodology 
and second quantization approach. Thus, the diagram 
should take into account the potential screening (includ¬ 
ing the polarization function frequency dispersion) and 
corrections to the electron propagator (self-energy). In 
this appendix we argue that the effects of interaction po¬ 
tential screening and the electron free energy are of the 
same order as the next perturbation theory correction. 

First, we neglect corrections to the electron prop¬ 
agator. We denote the free electron propagator as 
Go(co’, qi, q 2 ,i) {j is the LL in which the electron propa- 
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gates; qi /2 are the 2 -momenta — the electron propagator 
in uniform magnetic field is not translation-invariant, so 
it contains two momenta coming from the Fourier trans¬ 
forms with respect to two spatial coordinates). We de¬ 
note the photon propagator as D{uj, q) (it is translation- 
invariant, so depends only on one spatial momentum). 
Then the propagators can be expressed as follows (in¬ 
finitely small imaginary parts of the denominators are 
not written out explicitly): 


Go(w,qi,q2,i) 


/i(qi.q2) 

u-Ej 


(Cl) 


D(w,q) 


^>(q) 

1 + ^2i;(q)n(q, w) ’ 


(C 2 ) 


/j(qi,q 2 ) are smooth functions without poles, n(q,cj) 
is the polarization function, u(q) is the unscreened 
Coulomb interaction potential, the product u(q)n(q, cj) 
as a function of q is a smooth function without poles. 

The diagram we are interested in can be expressed then 
via the propagators (we will need only the uj dependence 
for our analysis, so the dependence on the spatial mo¬ 
menta is not written explicitly): 


Diagram = / duj G{uj,ji) x 

J over momenta J 

D{^E\ — uj)G{—lu + El + E2^ j2)D[—uj -h Ei + E 2 — -S 4 ), 

(C3) 

where Ei are the energies of the incoming and outgoing 
electrons. In our case both are in the same LL j: Ei = 
E¥'^ ^ Ji and j 2 denote the levels to which the electrons 
hop, just like in sub-subsection [TTl.2c[ 

All the functions (propagators) participating in the 
expression are the time-ordered correlation functions at 
zero temperature and their poles are situated in such a 
way that one can do Wick rotation uj ^ iuj (the zero 
of energies is the partially filled LL E^^^ = 0). Then, 
concentrating on the denominators, we can write the di¬ 
agram via the integral over Euclidean frequency: 


Diagram (x / / duj- -^ x 

J over momenta J — oo '^ji 

l + Pv{...)U{...,-iuj)' 

Evidently, the denominators of the electron propagators 
determine the region of frequencies which give dominat¬ 
ing contribution to the integral. And this region is the 
interval centered at zero frequency with the width of the 
order of the energy distance to the closest higher LL. 
Consider the difference of the full D{uj,q) and the bare 


Do(cj,q) = u(q) photon propagators: 


Do(w,q) -D(w,q) 


/^^(q)^n(q,u)) 

1 + ^2i;(q)n(q, w) 


D{lo, q) X O 




V itvhujc 


(C5) 


Therefore, replacing D(cc;,q) with Do(co’,q) brings in a 
relative error proportional to the ratio between the typi¬ 
cal interaction energy and the distance to the closest LL. 
Since we neglect next order perturbation theory correc¬ 
tion, taking the screening into account would exceed the 
accuracy of our calculation. In principle, it is possible 
to justify use of the screened interaction by expanding 
in other parameters (e.g. 1/N expansion [38]). However, 
in that case one needs to use not just the zero-frequency 
polarization function as in ( [ 6 ^ , but to take into account 
the polarization function frequency dispersion. We do 
not do this and use the unscreened interaction for the 
calculation. 

Now we discuss the effect of corrections to the electron 
propagator on the diagram. These corrections should 
lead to LLs shifting and broadening. It is natural to 
expect that the shifts and level widths are of the order of 
the Coulomb interaction energy scale /{nl). This would 
change the region of important frequencies by something 
proportional to the /{nl). 

Therefore, neglecting the effects of the polarization 
function and the electron free energy in the diagram in¬ 
troduces a relative error proportional to the ratio be¬ 
tween the typical interaction energy and the distance to 
the closest LL. Since we neglect the next order perturba¬ 
tion theory correction (and, moreover, the second order 
three-body interaction), taking these effects into account 
would exceed the accuracy of our calculation. 


Appendix D: Additional data on the Moore-Read 
PfafRan state stability 

In this appendix we provide additional details regard¬ 
ing the stability of the Moore-Read state. 

Eiguresj^ and[^ show the dependence of the overlap 
of the exact ground state of the system with the Moore- 
Read Pfafhan for A" = 12 particles on the magnetic field 
and the dielectric constant at 1/ = 50 meV for the (1,-1) 
LL and at U = 30 meV for the (2,-1) LL respectively 
with no virtual hopping to the nearby LLs taken into ae- 
eount. Screening is still taken into account. We do not 
show the data in the region where the perturbative anal¬ 
ysis is not applicable according to type S estimate (with 
typical interaction energy scale needed for the estimate 
taken to be the sereened potential zeroth pseudopotential 
at the LL under consideration). The region of inappli¬ 
cability of our theory according to the type C estimate 
(with typical interaction energy scale taken to be the bare 
Coulomb potential zeroth pseudopotential) is hatched. 
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Figure 6: Color plot of the overlap of the ground state 
with the Moore-Read PfafRan for 12 particles as a 
function of the magnetic field B and the dielectric 
constant k, with no virtual hopping corrections taken 
into account, (a) - for the (1,-1) LL at = 50 meV, 
(b) - for the (2,-1) LL at f/ = 30 meV. Contours show 
the lines of constant overlap. The region where perturbative 
analysis is not applicable according to the type C estimate is 
hatched. Data is not shown beyond the region where pertur¬ 
bative analysis is applicable according to the type S estimate. 
(Color online). 


As one can see (compare with Fig. [^, the effect vir¬ 
tual hopping corrections, defined in sub-subsection [lIL2c[ 
have on the phase diagram within the region of appli¬ 
cability of perturbative treatment of LL mixing is very 
significant, especially for the (1,-1) LL. 

Figures and[^ show the numerically found gap be¬ 
tween the exact ground state of the system with 12 par¬ 
ticles and its exact first excited state. In these figures, 
the gap is plotted as a function of the magnetic field and 
the dielectric constant at 1/ = 50 meV for the (1,-1) LL 
(Fig.[^) and at 1/ = 30 meV for the (2, —1) LL (Fig.[^) 
respectively. Virtual hopping to the nearby LL is taken 
into account. We do not show the data in the region 
where the perturbative analysis is not applicable accord¬ 
ing to type S estimate. The region of inapplicability of 
our theory according to the type C estimate is hatched. 
The blue lines are the overlap level lines from Fig. i Fig- 
ures and show the same data on the gaps in the 
units of typical Coulomb energy ^I[kI\ where I is the 
magnetic length I = ^Jhcj(eB). 

The behaviour of the gaps for the two levels is quite dif¬ 
ferent. For both levels the gaps generally increase as the 
magnetic field increases and as the dielectric constant de¬ 
creases. This can be partially (for magnetic field) or fully 
(for the dielectric constant) attributed to the increase of 
c^IkI. However, there are some interesting features in 
the plots. For the (1,-1) LL the gap drops to almost 
zero at a line in the hatched region. For the (2,-1) LL 
there are two such features: one in the hatched region 
and another slightly above the hatched region (better 
seen in Fig. than in Fig. [^). Although two of these 
features are situated in the hatched region where the ap¬ 
plicability of our methodology can be questioned, and 
the third one is much less pronounced, these feature, ob¬ 
tained for a finite number of particles, may suggest that 
a transition between different topological orders happens 
across those lines. Such transitions are likely to happen 
at small dielectric constants where the system physics 


Figure 7: Color plot of the gap between the ground 
state and the first excited state computed for 12 par¬ 
ticles as a function of the magnetic field B and the di¬ 
electric constant k,. (a) - for the (1,-1) LL at 1/ = 50 meV, 
(5) - for the (2,-1) LL at 1/ = 30 meV. The region where 
perturbative analysis is not applicable according to the type 
C estimate is hatched. Data is not shown beyond the re¬ 
gion where perturbative analysis is applicable according to 
the type S estimate. The blue lines coincide with the overlap 
level lines from Fig.(Color online). 



Figure 8: Color plot of the gap between the ground 
state and the first excited state computed for 12 par¬ 
ticles as a function of the magnetic field B and the di¬ 
electric constant k. (a) - for the (1,-1) LL at 1/ = 50 meV, 
(5) - for the (2,-1) LL at 1/ = 30 meV. The region where 
perturbative analysis is not applicable according to the type 
C estimate is hatched. Data is not shown beyond the re¬ 
gion where perturbative analysis is applicable according to 
the type S estimate. The blue lines coincide with the overlap 
level lines from Fig.(Color online). 


becomes essentially non-single Landau level. 

We remind the reader that for the (2,-1) LL the high 
overlap with the Moore-Read Pfaffian state is achieved 
near H = 8T at > 20, as can be seen from Fig. 

Figures!^ and[^ show the dependence of the gap on 
the magnetic field at 7/ = 30 meV and /^ = 20 and 25 for 
the (2, —1) LL in K and in /{nl) units respectively. 

Figures andpRp show the dependence of the gap on 
the dielectric constant at 7/ = 30 meV and fixed H = 8 T 
for the (2, —1) LL in K and in /{nl) units respectively. 

Now we give some data on overlaps and gaps from nu¬ 
merical diagonalization for different numbers of particles 
N = 8,10,12,14. 

Figure [TT1 compares the data on the overlap with the 
Moore-Read Pfaffian and the gap to the first excited state 
for the (1,-1) and the (2,-1) LLs in BLG and the non¬ 
relativist ic n = 1 LL. The data are shown for different 
numbers of particles N = 8,10,12,14. For the (1,-1) 
and the (2,-1) LLs in BLG external parameters are set 
to be near the maximum overlap regions (according to 
Fig-lU) at several values of the dielectric constant k] we 




















24 



Figure 9: Dependence of the gap between the ground 
state and the first excited state computed for 12 par¬ 
ticles at the (2,-1) LL for U = 30 meV and k, = 20 
(blue dashed line) and 25 (black solid line) as a func¬ 
tion of the magnetic field B. (a) - in K, (b) - in e‘^/( kI) 
units. Only the part where perturbative analysis is applicable 
according to the type S estimate is shown. (Color online). 



Figure 10: Dependence of the gap between the ground 
state and the first excited state computed for 12 par¬ 
ticles at the (2,-1) LL for U = 30 meV and B = 8 T as 
a function of the dielectric constant k. (a) - in K, (b) 

- in /{nil) units. Only the part where perturbative analysis 
is applicable according to the type S estimate is shown. The 
cusps present in the plots are due to the fact that simulations 
were performed at a discrete set of dielectric constants Hi with 
a step 6Hi = 1.5. Therefore, the data accurately represents the 
region where the dependence is smooth, while where the gap 
changes quickly the true data between the cusp points may 
deviate from the lines drawn. 


also present the data for 5 = 8, Hi = 2b for the (1,-1) 
LL. The gap is presented in units of the typical Coulomb 
energy /{nl). The overlaps for the (2, — 1) LL, as well as 
for the (1, — 1) LL at /^ = 40, appear to be more stable as 
the number of particles grows than for the non-relativistic 
system. It is interesting to note that for the (2,-1) LL 
the overlap is slightly higher at = 12 particles than 
at A/" = 8,10, and 14, while for the (1, —1) LL it is vice 


Figures present the dependence of the overlap 

and the gap at several parameter points for the (2,-1) 
LL. Points are taken to be at the same mini-gap and the 
dielectric constant values as in the previous figure, but 
the magnetic field changes. We took the points in the 
maximum overlap region, slightly to the left and to the 
right of it. 


Figure 11: Dependence of the overlap with the PfafRan 
and gap to the first excited state on the number of 
particles N (N = 8,10,12,14). Solid black line is for the 
non-relativistic n = 1 LL. Dashed lines are for the (2, —1) LL 
at = 30 meV: ai k = 2b and B = 8 T (black), at = 20 
and B = 7.5 T (blue). Dot-dashed lines are for the (1,-1) LL 
at P = 50 meV: at = 25 and B = 8 T (black), at = 40 
and B = ST (green). The gap is presented in units of typical 
Coulomb energy e^/(/^/). (Color online). 
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Figure 12: Dependence of the overlap with the Moore- 
Read PfafRan and gap to the Rrst excited state on the 
number of particles N {N = 8,10,12,14) for the (2, —1) 
LL at 1/ = 30 meV, k, = 20. Shown are the dependences for 
B = 7.b T (dashed black line), B = 6 T (dashed grey line), 
and B = 9 T (dashed purple line). (Color online). 
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Figure 13: Dependence of the overlap with the Moore- 
Read PfafRan and gap to the Rrst excited state on the 
number of particles N {N = 8,10,12,14) for the (2, —1) 
LL at P = 30 meV, k, = 25. Shown are the dependences for 
B = 8 T (dashed black line), B = 6.b T (dashed grey line), 
and B = 9.b T (dashed purple line). (Color online). 
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